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Abstract. We present high resolution spherically symmetric relativistic magnetohydrodynamical simulations of 
the evolution of a pulsar wind nebula inside the free expanding ejecta of the supernova progenitor. The evolution 
is followed starting from a few years after the supernova explosion and up to an age of the remnant of 1500 years. 
We consider different values of the pulsar wind magnetization parameter and also different braking indices for the 
spin-down process. We compare the numerical results with those derived through an approximate semi-analytical 
approach that allows us to trace the time evolution of the positions of both the pulsar wind termination shock and 
the contact discontinuity between the nebula and the supernova ejecta. We also discuss, whenever a comparison is 
possible, to what extent our numerical results agree with former self-similar models, and how these models could 
be adapted to take into account the temporal evolution of the system. The inferred magnetization of the pulsar 
wind could be an order of magnitude lower than that derived from time independent analytic models. 
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1. Introduction 

Pulsars are rapidly spinning magnetized neutron stars 
that usually form as the result of the core collapse of 
massive stars (8 — 16 Mq) in supernova events (SN). 
The typical energy released in a supernova explosion is 
of order ~ 10 53 erg. Most of this is carried away by 
neutrinos, while only a small fraction (about 1%) goes 
into a blast wave that sweeps up the outer layers of the 
star and produces a strong shock propagating in the sur- 
rounding medium. The ejected material is initially heated 
by the blast wave and set into motion. Then, while the 
heat is converted into kinetic energy, the ejecta acceler- 
ate until the pressure becomes so low as to be dynami- 
cally unimportant. When this happens the material finally 
sets into homologous expansion IjChevalier fc Soker~1 989 
IMatzner fc McKee 1999(1 . This phase is usually referred to 
as free expansion of the ejecta. 

As a consequence of the electromagnetic torques acting 
on it, the pulsar supplies a late energy input to the rem- 
nant in the form of a relativistic magnetized wind, mainly 
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made of electron-positron pairs and a toroidal magnetic 
field l|Colrireich fc Julian 19691 IMichel fc LiT999ll . Most 
of the pulsar rotational energy is carried away by this 
wind, whose propagation velocity is ultra-relativistic, with 
typical Lorentz factors that, far enough from the light 
cylinder, are estimated to be in the range 10 4 - 10 7 . The 
interaction of the wind with the ejecta expanding at non- 
relativistic speed produces a reverse shock that propa- 
gates toward the pulsar (Re es fc Gunn 1974jl . In the re- 
gion bound by the wind termination shock on the inner 
side, and by the ejecta on the outer side, a bubble of rela- 
tivistically hot magnetized plasma is created. This shines 
through synchrotron and Inverse Compton emission in a 
very broad range of frequencies, from radio wavelengths 
up to gamma rays: this is what we call a pulsar wind neb- 
ula (PWN) or plerion. 

The evolution of a PWN inside the free expanding 
ejecta depends on many different parameters such as 
the pulsar luminosity, the density and velocity distri- 
bution in the SN ejecta IjDwarkadas fc Chevalier 1 998 
IFeatherstone et al. 20011 IBlondin et al. 1996)l . the pres- 
ence of large and/or small scale inhomogeneities 
IIChevaHer fc Soker 19891 |Campbell et al. 2003| ). In the 
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case of constant luminosity, and if spherical symmetry 
is assumed, it is possible to derive a simple evolutionary 
equation for the radius of the PWN as a function of time, 
that results in a power law IjChevalier fc Fransson 1992 
Ivan der Swaluw et al. 200l|l . In the case of SN ejecta with 
a constant density profile the PWN contact discontinuity 
evolves as i 6 / 5 . For a more detailed description of the var- 
ious phases of the PWN-SNR evolution see, for example, 
van der Swaluw et al. (2001) and references therein. 

While many analytic and numerical models ex- 
ist in the literature for the evolution of SNRs, un- 
til recently only two classes of analytic models in the 
proper relativistic magnetohydrodynamical regime have 
been presented for PWNe: the steady state solution by 
Kennel fc Coroniti (1984)| (KC hereafter), and the self- 



similar solution by |Emmering fc Chevalier (1987) (EC 
hereafter), which allows for a non zero velocity of the ter- 
mination shock. 

Only lately the evolution of PWN-SNR systems 
has began being investigated through numerical sim- 
ulations. These have been performed mainly in the 
classical hydrodynamical (HD) (Blon din et al. 2 001 
Ivan der Swaluw et al. 200Tjl . or classical MHD 
Ijvan der Swaluw 2003(1 regime. However the recent 
development of codes for relativistic magnetohydrody- 
namics (RMHD) allows one to investigate such systems 
in a proper regime and to quantify the accuracy of 
approximate analytic solutions (Buccian tini et al. 2003(1 . 

Both KC and EC models rely on two strong assump- 
tions: a constant pulsar spin-down luminosity, and a con- 
stant velocity at the outer boundary of the nebula, nei- 
ther of which applies to a real case nor is consistent 
with the PWN evolution inside free expanding ejecta. 
While both assumptions are known to be unrealistic, the 
most crucial one, as far as the long-term evolution of 
the system is concerned, is probably that of constant 
pulsar energy input. As we have already mentioned, the 
PWN is powered by the rotational energy lost by the star 
due to electromagnetic braking, and this loss translates 
into an increase with time of the pulsar rotation period 
flLyne fe Graham-Smith 1998Q . 

In the case of a dipolar magnetic field the torque 
exerted on the star results in the following relation be- 
tween the spin-down rate and the pulsar frequency f2 (e.g. 
IMichel fc Li"T999)l : 



(1) 



while the power supplied to the wind changes with time t 
as: 

La 



IQ.Q. = L(t) = 



(2) 



where / is the momentum of inertia of the pulsar, r is 
the characteristic spin-down time, and L Q is the initial 
pulsar luminosity. More generally, if the field is not exactly 
dipolar one can write the pulsar (or wind) luminosity as: 

Lit) = - — ^r— , (3) 
w (1 + t/r) n V ' 



where n=(/3 + l)/(/3— 1), with (3 the braking index. 

Estimated values of L Q may be up to 10 38 -10 40 erg/s. 
Determining the braking index from observations is 
extremely complicated, as it requires detailed and 
precise pulsar timing over long time-spans. A mea- 
sure of n is presently available for four pulsars only 
l(Camilo et al. 2 000). Among these, one, the Vela pulsar, 
has n = 6, while all the others have 2 < n < 3. The value 
of t can be derived from the pulsar period and spin-down 
rate once n is known. 

In the following we present ID high resolution numer- 
ical simulations of the structure and evolution of a PWN 
inside free expanding SN ejecta. We extend, for this ini- 
tial stage, that lasts for about 2000-3000 years (until the 
reverse shock propagating in the SNR collapses on the 
PWN), previous work on the subject by Bucciantini et al. 
(2003). 

The plan of the paper is as follows. In Section [21 we 
present a semi-analytic model for the evolution of the 
PWN radius and the pulsar wind termination shock, valid 
for any power-law profile of the ejecta and for any value 
of n. After briefly describing, in Section [3 the numerical 
code employed and the initial conditions of our simula- 
tions, we present, in Section the results obtained for 
different values of the braking index, n = 0,2,3, and dif- 
ferent magnetizations of the wind. We discuss these results 
and compare them with the expectations of the analytic 
EC model. A formula is derived for the spin-down factor- 
ization, and the new results are applied to the case of Crab 
Nebula. In Section El we summarize our conclusions. 



2. A semi-analytic preamble 

In this section we shall derive some approximate relations 
for the evolution of a PWN interacting with the SN ejecta: 
these will be of use to interpret the numerical results of 
Section^ Our main goal will be to find an expression for 
the time evolution of the position of the termination shock 
(Rts(t)) and of the contact discontinuity (R c d(t))- 

First of all, using the results found by Bucciantini et al. 
(2003), we write the total (magnetic plus particle) energy 
inside the PWN as: 



(4) 



with P c d the total pressure at the contact discontinuity. 
The validity of this equation is easy to prove in the two ex- 
treme cases: if we assume for the PWN a non magnetized 
bubble with constant thermal pressure or a magnetically 
dominated bubble with magnetic field B oc r . However, 
the simulations we present in the following (Section [3} 
support the idea that this relation holds to a very good 
approximation (within 0.3% error), whatever the spatial 
distribution of magnetic field and thermal energy in the 
nebula, and for all values of n considered. 
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From Eq. it follows that the evolution of E pwn , in- 
cluding the pulsar input and the effects of adiabatic ex- 
pansion, can be written in the form: 



^(4tt Rl d P cd )=L{t) 



4tt R z cd P cd R cd . 



Integration by parts leads to: 



Pcd(t) = 



1 



ATrR cd (ty 



L(s)R cd (s)ds 



(5) 



(6) 



This latter equation makes it clear that, for given ejecta 
properties, once L(t) is assigned, both R c d{t) and P c d(t) 
are uniquely determined, independently of the wind mag- 
netization. 

What does depend on the wind magnetization is the 
internal pressure profile: while in a HD case the pressure 
just behind the shock would be almost equal to P c d(t), 
in a MHD case it will be higher. Moreover the ratio of 
the post-shock pressure over the wind ram pressure (just 
upstream of the shock) will depend on the termination 
shock speed. Approximate pressure balance at the shock 
will still hold in all cases where the shock speed is low 
compared to the speed of light. We can then write: 



L(t) 



^Rl(t)c 



Pu(t) 



Ptsjt) Jo L(s)R cd (s)ds 
P cd {t) 4vri?4 d (t) 



(7) 



where in the latter equality we made use of the expression 
in Eq.|H|for P cd (t). It should be noticed that both P cd and 
Pts arc here defined as including only the contribution of 
thermal plus magnetic pressure in the nebula, whereas the 
particle ram pressure is neglected. The latter is expected 
to be negligible at the outer boundary of the nebula, while 
in a truly steady-state situation it would contribute 1/3 
of the total pressure just behind the shock. 

Eq. Hallows one to derive the time-evolution of the ter- 
mination shock radius, once R c d(t) is known, and a model 
for the behavior of the ratio Pts /Pcd as a function of time 
is provided. As far as the first task is concerned, a sim- 
ple analytic description of how the contact discontinuity 
position evolves with time can be found under the thin- 
layer approximation. In this approximation, R c d{t) can be 
derived from the following equation: 



R cd (t')L(t') dt' = R 2 cd [M sh R cd + M sh [R cd 



to express the latter in terms of the enclosed mass, in the 
following way: 



M 



W 



(a-l)(a+l)(a + (3-0(a-l)) 



3-? 



(9) 



We have chosen the factor W in the above equation in 
such a way that the following simple expansion law holds 
for a constant energy input by the pulsar: 



(10) 



where a = (6 + £)/(5 + £) , and the latter equality defines 
R cd . In the general case of a fading L{t) (described by 
Eq. |3J) , the evolution of i? c( j is not exactly a power-law. 
While, as shown by Eq. ^3 at times much smaller than r 
the best- fit power-law expansion has an index a = a , at 
very late times the evolution of R c d is well approximated 
by a linear law (i.e. a = 1). 

We employ a series expansion that is capable of de- 
scribing the evolution of R cd at all times, and which has 
a functional form that allows the time- integral in Eq.Qto 
be performed analytically. Let us introduce the variable 
s = (t/r)/(l + t/r): the temporal range [0, oo] maps into 
the range [0, 1] for s. 

Let us then express R c d{t) by the series: 



ctl 



— Y 

j=0 



c 4 s 



(11) 



with Co = 1. The factor in front of the series guarantees 
that, for iCr, Rcdifi} is well approximated by R c d,o(t), 
while at larger times R c d(t) oc t. 

An advantage of this functional form is that the series 
converges at all times, while a simple power series of t has 
been found not to converge for t > r. The values of the 
various coefficients can be easily obtained in the case of 
constant L (i.e. n = 0), where R c d follows Eq. which 
implies: 



i=0 



(1 



E 



(-l)*r(2-a ) 



^Vr(2-a 



i)T(l + i) 



(we have used the binomial expansion). 

In the case of a general n, Eq.[S]must be directly solved 
in order to obtain the coefficients c;. The first ones are: 



^i)j.(S) 



where M s h is the mass of the shocked 
ejecta. Eq. |S] is obtained by combining (e.g. 
|Reynolds &l Chevalier 1984| )mass and momentum conser- 
vation with the energy conservation law of Eq. [S] 

It is possible to derive a power series approximation of 
Eq.|Sl which allows an analytic description of the evolution 
of i? c( j at any time. It is worth stressing that the series 
expansion we present in the following can be applied to 
the case of a general braking law for the plerion-feeding 
pulsar (Eq. I3J) , as well as to the case of a general power- 
law density profile of the ejecta p e j °c t^~ 3 r~^. We prefer 



ci 



C-2 



(49-9£)-n(ll-20 



(12) 



245 - 94£ + 9£ 2 
- ((49 - 9O 2 (-7038 + 4029£ - 764£ 2 + 48£ 3 )) 

-n(-6398469 + 5961745£ - 2216513£ 2 + 411061£ 3 
-38028£ 4 + 1404£ 5 ) + n 2 (624393 - 577770£ 

-213273C 2 - 39259£ 3 4- 3604£ 4 - 132£ ! 
(2(49 - 90 2 (5 -£) 2 (H< 7 3 - 476£ + 48£ 2 ) , (13) 



while the further ones are too complicated to be listed 
here. However, it can be seen that the coefficient c, is a 
polynomial of i-th degree in n, and in Table ^ f° r a few 
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Table 1. The first few coefficients Cj for the series expansion in Eq. El The polynomial expression of a as a function 
of the pulsar spin-down index is given for different power-law profiles of the ejecta. 





= 


Cl 




0.200000 


- 0.044898 n 










C2 


= 


0.120000 


- 0.045438 n + 0.004434 n 2 










C3 


= 


0.088000 


- 0.043397 n + 0.007204 n 2 


- 0.000383 n 3 








C4 


= 


0.070400 


- 0.041023 n + 0.008888 n 2 


- 0.000802 n 3 


+ 0.000024 n 4 




= i 


Cl 




0.250000 


- 0.056250 n 










C2 


= 


0.156250 


- 0.059742 n + 0.005875 n 2 










C3 


= 


0.117188 


- 0.058696 n + 0.009905 n 2 


- 0.000540 n 3 








C4 


= 


0.095215 


- 0.056620 n + 0.012563 n 2 


-0.001176 n 3 


+ 0.000037 n 4 




= 2 


Cl 




0.333333 


- 0.075269 n 










C2 




0.222222 


- 0.086234 n + 0.008582 n 2 










C3 




0.172840 


- 0.088634 n + 0.015326 n 2 


- 0.000864 n 3 








C4 




0.144033 


- 0.088316 n + 0.020295 n 2 


- 0.002000 n 3 


+ 0.000069 n 4 




= 3 


Cl 




0.500000 


- 0.113636 n 










C2 




0.375000 


- 0.149268 n + 0.015168 n 2 










C3 




0.312500 


- 0.166823 n + 0.030039 n 2 


- 0.001785 n 3 








C4 




0.273438 


- 0.176633 rc + 0.042976 n 2 


- 0.004568 n 3 


+ 0.000177 n 4 



choices of £, we list the values of these coefficients up to 
the 4 th order. 

We have tested these approximate analytic solutions 
with a numerical model, in thin-layer approximation, for 
the case of flat ejecta (£ = 0), which is the one relevant for 
the situations considered here. The discrepancy increases 
with time up to an asymptotic value of 18%, 13%, 11% 
and 9%, respectively, for the first to fourth degree ap- 
proximations. In the range of times shorter than 3r (as in 
our simulations, see next section), the errors are not larger 
than 4.7%, 2.2%, 1.2%, and 0.7%, for approximations of 
increasing degree. 

The expression chosen for R c a(t) fEa. lllfl contains only 
terms proportional to (i/r) M (l + t/r) v , and therefore the 
integral J LR dt can be evaluated as a series of hypergeo- 
metric functions: 



L(t')R(t')dt' = R cd L T 



^ 1 

i=0 



t 

l \ T 



l+a +i 



2-Fi (n - 1 + a Q + i, 1 + a Q + i, 2 + a Q + i; —t/r) , (14) 

where 2,Fx(a, b, c; z) is the hypergeometric function. 

Given Eq. ^2 Eq. El an( i the coefficients Cj , what is 
left to find, in order to obtain the time evolution of Rt s 
from Eq. is an expression for the ratio Pts/Pcd as a 
function of time. This can be easily accomplished under 
the approximations that the PWN always adjusts itself to 
a quasi-steady solution and that the fluid motion is non- 
relativistic everywhere behind the termination shock. Let 
us define the quantity II = P/(p u c 2 j 2 ), where p u is the 
matter density and 7 the Lorentz factor of the realtivistic 
wind at the termination shock. Following KC, one finds 
for II the approximate expression: 



1%) 



27 



G(y) 4 



3y 2 



(15) 



with y a non-dimensional coordinate related to the dis- 
tance from the shock r, and to the magnetization of the 
wind a (supposed to be low), 



y(r) 



81 a r 



2 Rts 

and the function G(y) (see KC) given by: 



(16) 



G(y) = 1 + [l + y 2 + Va + y 2 ) 2 -! 

1/3 



-1/3 



+ l + y 2 + v /(l + y2)2_ 1 



(17) 



We use for the magnetization parameter a the definition 
a = B 2 / {Airpc 2 ^ 2 ) with B and p the wind magnetic field 
and matter proper density respectively. 

We emphasize that the normalized pressure profile in 
Eq. E| is straightforward to obtain as a solution of the 
steady-state MHD equations under the assumption that 
the post-shock bulk Lorentz factor of the fluid is 7 = 1 
and that a <C 1. The first term in Eq. El is the thermal 
pressure, while the second term is the magnetic contribu- 
tion, as can be easily checked evaluating the expression in 
the proper limits of a and hence y. 

The ratio Pts/Pcd in Eq. [7| which we shall indicate 
hereafter as 1/K, can be expressed in terms of II as: 



1 



Pts(t) 



n„ 



K(t) P cd (t) n(y(R cd (t))) 



(18) 



with n = U(y(R ts (t)) = 11(^81(7/2). It is apparent that 
under these simplified assumptions the ratio between the 
value of the pressure at the termination shock and that at 
the contact discontinuity depends on time only implicitly, 
through the ratio between Rt s and R c d- 
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In light of Eq. ^3] we can rewrite Eq. [7| as an implicit 
equation for y cd {t) = y(R cd (t)) = y/%\<j/2 R cd {t) / R ts {t): 

y 2 cd (t) U(y cd (t)) « ^-11 H Q(t) , (19) 

where Q(t) can be determined using Ea. II II and Ea. 1141 



Q(t) 



Cl (t/r) 



^(l + «„ + i) 
F\ [n - 1 + a a + i, 1 + a Q + i, 2 + a Q + i; 



(20) 



The right-hand side of Eq. ^|is then a known function of 
time alone. Hence the equation can be easily solved numer- 
ically to obtain the shock position as a function of time. 
As we shall see in the following, this approach, despite be- 
ing extremely simplified, allows one to trace the evolution 
of the termination shock radius with an accuracy of order 
15%, once the value of R cd and R ts are known at a refer- 
ence time. In fact, as already mentioned, the first equality 
in Eq.[7]holds only within a 15-25%, and this means that 
Eq. 1191 must be normalized to observations (or to simu- 
lations) to provide a good estimate of the position of the 
termination shock at any given time. 

3. Numerical simulations 

The simulations we present have been performed by us- 
ing the newly developed scheme by Del Zanna et al. 
IjDel Zanna fc Bucciantini 20021 IDel Zanna et al. 2003|l . 
We refer the reader to the cited papers for a detailed de- 
scription of the code, and of the equations and algorithms 
employed. This is a high resolution conservative (shock- 
capturing) code for 3D-RMHD based on third order ac- 
curate ENO-type reconstruction algorithms. The approx- 
imate Riemann solver employed is the two-speed HLL 
flux formula, which does not make use of time-consuming 
characteristics waves decomposition. Given the presence 
of very strong shocks we have used second order recon- 
struction, to reduce post-shock oscillations. 

We have used a single fluid model, assuming an adia- 
batic coefficient equal to 4/3 also for the SN ejecta. This 
makes the ejecta more compressible but does not change 
the temporal evolution of the PWN. We have chosen to use 
a single fluid because, in numerical simulations, the use of 
two different adiabatic coefficients on a contact disconti- 
nuity with a very large density jump (density may change 
by factors of order 10 6 — 10 7 ), leads to the formation of 
spurious disturbances that tend to propagate back into 
the PWN flShyue 1998| IKarni 19981 IKun fc Jishan 19981 
IBucciantini et al 2 003 ) . 

When a magnetic field is present the speed of a shocked 
wind cannot drop to zero but tends to a finite value V aS y 
(see e.g. KC). If the contact discontinuity velocity is close 



to this asymptotic value, even small fluctuations, originat- 
ing at the PWN-SNR interface, can produce substantial 
variations in the internal structure. 

Another source of problems is numerical diffusion at 
the contact discontinuity itself. This has the effect of 
spreading the density jump at the contact discontinuity 
over a few computational cells. A criterion needs then to 
be established for the identification of the position of the 
contact discontinuity (R cd ), in order to compare the simu- 
lation results with the existing analytic models. We found 
that a convenient choice is to identify R cd with the posi- 
tion where the fluid velocity is equal to the contact discon- 
tinuity velocity: v(R cd ) — V cd . This offsets R cd by 3-4 % 
(R cd is 3-4 simulation cells further out) with respect to 
the radius at which the density jump begins. 

3.1. Initial conditions 

Simulations have been performed on a logarithmic radial 
grid, with 200 cells per decade. This allows one to resolve 
with sufficient accuracy the inner region so that the ter- 
mination shock always remains inside the computational 
domain and injection ambiguities such as those discussed 
by Bucciantini et al. (2003) are avoided. At the same time, 
with our choice of the grid we are able to follow the sys- 
tem evolution from very early times (a few years) after the 
SN explosion and up to an advanced age, maintaining dif- 
fusion effects homogeneous throughout the evolution. We 
set continuous conditions at the outer boundary (zeroth 
order extrapolation). No radiation cooling is included. 

For the free expanding ejecta we have chosen the fol- 
lowing profiles (Chev alier fc Fransson 19 92): 



Pcj = At 
v = r/t = V r/R 0l 



R Q {t) = V t, 



(21) 



with A = 8.7 x 10 6 g s 3 /cm 3 and V Q = 5.27 x 10 8 cm/s, 
corresponding to an energy release in the SN explosion 
E = 10 51 erg and ejecta having mass M ~ 4M Q . The 
pulsar wind is created injecting mass, momentum-energy, 
and a purely toroidal magnetic field in the first computa- 
tional cell, with a total luminosity that depends on time as 
described by Eq. [31 with L — 5 x 10 39 erg/s, and r = 500 
years. Three different values for the spin-down index have 
been used: n = 0, 2, 3. No magnetic field is initially present 
in the SNR nor in the ISM. 

The simulations are initialized with a 5 years old PWN 
surrounded by a thin shell of swept up ejecta at a radius 
Red = 0.022 ly. This shell contains the ejecta material re- 
moved from the origin by the relativistic bubble. The evo- 
lution of the PWN is followed up to an age t = 3 r = 1500 
years. The internal profile of the fluid in the PWN is taken 
from the EC solution with proper magnetization, fitted to 
the contact discontinuity velocity. After a short transient 
phase of about 10-15 years, partly due to numerical ef- 
fects, the nebula relaxes to a stable configuration. Our 
choice of initial conditions is such that the system is not 
far from the self-similar solution. This allows us to avoid 
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the long term transient that is otherwise observed when 
the pulsar wind is not switched on at the same time when 
the SN goes off IjChevalier fe Fransson 19921 Uun 1998|l . 
Self-similar solutions only exist when a is low enough so 
that the asymptotic velocity of the shocked wind is smaller 
than the contact discontinuity initial speed, V asy < V cc i- 
Three wind cases have been simulated: 

— Purely hydrodynamic wind (a = 0); 

— Weakly magnetized wind (a — 0.0016); 

— "Highly" magnetized wind (a = 0.003), with highly 
meaning close to the maximum compatible with the 
existence of a self-similar solution. 

In all cases the wind has a Lorentz factor 7 = 100, and 
p/pc 2 = 0.01. 

4. Numerical results 

4.1. Constant luminosity case 

The first step we take is the comparison of our numer- 
ical results with the analytic model by EC in the case 
with n = 0. The EC model is completely determined once 
the wind quantities (density, pressure, magnetic field and 
Lorentz factor) and the termination shock velocity are 
known. Given the value of the termination shock speed, 
Vtsi the shock jump conditions are evaluated at the termi- 
nation shock radius, Rt s , and then the self-similar, spher- 
ically symmetric, RMHD equations are integrated using 
the post-shock values as initial conditions. A singularity 
appears in the solution at some distance from the termi- 
nation shock: the position of this singularity corresponds 
to the outer boundary of the nebula. 

The self-similar solution of EC requires Vt s to be con- 
stant and equal to a given fraction of the contact discon- 
tinuity velocity V c d, also constant: 

//,ST - (22) 



CD and TS evolution 



v t 



R, 



-V cd . 



We notice that the latter equation can be derived di- 
rectly from the general expression in Eq. Let us assume 
a constant pulsar energy input, L = Lq and the evolution 
of R c d described by a power-law: R c d oc t a . The time- 
derivative of Eq. [3 gives: 



V ti 



Rt 
2R, 



ill 



3V cd 



Red RcdK 

~T K 



(23) 



which, in the self-similar case (i.e. a = 1 and dK/dt = 0), 
reduces to Eq.|221 More generally, for R c d(t) still described 
by a power-law but with an index a ^ 1 , if the variation of 
K in time can be neglected, we find R ts octO- 1 )/ 2 . While 
in a HD case, one will have K rs const s=s 1 (see also Ea.1151 
in the limit y — > 0), the time- variation of K will become 
more and more important the larger the magnetization of 
the wind. 

In Fig. ^ we show the evolution of R c d derived 
from our simulations: this turns out to be indepen- 
dent of the magnetization, as stressed by Bucciantini 




Fig. 1. Evolution with time for the n — case of the 
position of the contact discontinuity (solid line) and of the 
termination shock for the hydrodynamics (dotted line), 
the a = 0.0016 (dashed line) and the a = 0.003 (dash- 
dotted line) cases. 

et al. (2003). The temporal evolution agrees with the 
behavior R c d oc t 6 / 5 predicted by the analytic models 
HChevalier Fransson 1992llva,n der Swaluw et al. 2001 II . 
The value of Vts changes with the magnetization and is 
lower for larger values of a. 

We find that the evolution of the termination shock is 
well described by a power law in time Rt s oc t s , but the 
exponent changes with the wind magnetization: in the HD 
case we find S — 13/10, as predicted by Eq. [23] for K = 
and a = 6/5, while the value of 8 increases to 1.38 and 1.43 
in the cases with a — 0.0016 and a — 0.003 respectively. 
The evolution of the termination shock radius is very well 
described in these cases by Eq. ^5] for the proper value of 
the magnetization. 

We have verified that, when computing the appropriate 
EC model for comparison, using the value of V ts derived 
from the simulations, or that obtained from the solution 
of Eq. ^5] does not improve substantially the result with 
respect to using Eq.[53for given values of R c d, V c d and R ts - 
A general advantage that Eq.|221offers is that it allows one 
to estimate V* s from quantities that can all be measured 
directly, at least in principle. 

In Fig. [21 we compare the radial profiles of density 
and pressure (both thermal and magnetic) derived from 
the simulations with those computed based on EC. The 
agreement of the EC profiles with our results is extremely 
good and fails (see the pressure profile) only near the outer 
boundary of the nebula. This is a consequence of the self- 
similarity imposed in the EC model which leads to an 
unphysical singularity at the border. 

We want to point out some differences between the 
EC and KC models which lead us to conclude that the 
EC model is better suited as the basis for a comparison 
with the numerical results presented here. First of all, the 
EC model reproduces the positive pressure gradient in the 
post-shock region observed in the simulations (the pres- 
sure initially increases up to a value that is about 20-25 % 
larger than that immediately behind the shock), while 
in KC the pressure is always monotonically decreasing. 
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Density at t=1 500 yr 



CD and TS evolution 




Thermal pressure at t = 1 500 yr 




Magnetic pressure at t = 1 500 yr 



a) 0.100 r 



& 0.010 r 




R(ly) 



Fig. 2. Comparison of the results with EC (solid line) for 
the no spin-down case at t=1500 yr. Density, thermal pres- 
sure and magnetic pressure are shown for all the mag- 
netization we have considered: a = 0.003 (dotted line), 
a = 0.0016 (dashed line), a = (dash-dotted line). The 
total pressure at the border is the same in the various 
cases. We want to stress the the radius of the contact dis- 
continuity is about 3-4% greater that the position where 
the density jump starts. This explain differences with re- 
spect to Fig. 



Moreover EC gives a larger flow speed in the asymptotic 
region than KC, in better agreement with the simulation 
results. 

However, the comparison with the EC model is less 
satisfactory if we consider the size of the nebula. In com- 
paring the model with observations, the standard way to 
estimate the magnetization parameter is as follows: one 
finds the appropriate a so that the theoretical fluid speed 
at a distance from the shock corresponding to the ob- 
served value of R c d matches the boundary speed derived 



Fig. 3. Evolution with time for the n — 3 case of the 
position of the contact discontinuity (solid line) and of the 
termination shock for the hydrodynamics (dotted line), 
the a = 0.0016 (dashed line) and the a = 0.003 (dash- 
dotted line) cases. Now the evolution of the termination 
shock cannot be described as a power law in time. 

from observations. However, the values R(v = V c d) / Rt s for 
the EC model are lower than the values R c d/Rts in our 
simulations, and the discrepancy increases with a (Fig. 0) . 

This discrepancy is most likely a consequence of the 
changing boundary speed and is enhanced when the con- 
tact discontinuity moves with a velocity close to V asy . 
This is why in the HD case, when V asy = (and hence 
Vcd /V aS y 3> 1) the difference is very small. On the other 
hand, we see in Fig. that, if instead of matching the 
velocities, we consider the size of the nebula in the EC 
model as the radius at which the EC solution has a sin- 
gularity, we observe a much better agreement. Still some 
discrepancy remains in the magnetic cases but now it is 
well below 10 %. 



4.2. Cases with spin-down 

Let us now turn our attention to the main assumption 
underlying the EC solution, i.e. that of a constant pulsar 
luminosity. As previously noted iBucciantini et al. 20 03). 
the spin-down process has the effect of reducing the ram 
pressure in the wind and, as a consequence, we expect 
to find a ratio R c d/Rts greater than in the case with no 
spin-down. 

When the effects of the pulsar spin-down are included, 
the evolution of R c d(t) can no longer be described in terms 
of a fixed power law in time. Despite this, we find that 
the variation of the exponent is small enough that we can 
continue to approximate it as a constant whose value, as 
derived from the simulations, is found to be slightly less 
(~ 1.13 — 1.1) than in the case n = 0. A correct descrip- 
tion of the time-evolution of R c d can still be derived as 
discussed in Section [21 Eq- IHI with the appropriate val- 
ues of the coefficients Cj, as reported in Tabled is found 
to provide a very good approximation for R c d{t) within 
the uncertainties discussed in the same session. The same 
is true for Rt s (t): this can be computed as described in 
Section with an accuracy of order 10-15 %, with the 



8 



N. Bucciantini et al.: The effects of spin-down on the structure and evoiution of PWNe 



Density at t=1 500 yr 
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Magnetic pressure at t = 1 500 yr 
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Fig. 4. Comparison of the results with EC (solid line) for 
the n — 2 case at t=1500 yr. Density , thermal pres- 
sure and magnetic pressure are shown for all the mag- 
netizations we have considered: a = 0.003 (dotted line), 
cr = 0.0016 (dashed line), a = (dash-dotted line). The 
total pressure at the border is the same in the various 
cases. Comments are the same as for Fig. |21 



Fig. 5. Comparison of the results with EC (solid line)for 
the n — 3 case at t=1500 yr. Density, thermal pressure and 
magnetic pressure are shown for all the magnetizations 
we have considered: a = 0.003 (dotted line), a = 0.0016 
(dashed line), a = (dash-dotted line). The total pressure 
at the border is the same in the various cases. Comments 
are the same as for Fig. |21 



error increasing with increasing n and a. The reason for 
this can be easily understood: our approach for the com- 
putation of Rt s {t) is based on the assumption that the 
evolution of the nebula proceeds slowly. This becomes an 
increasingly bad approximation as n and a increase. 

In Figs.0]and[S]we plot the results of our simulations in 
the case n — 2 and n — 3 for all the different values of a we 
adopted. A comparison is made between the profiles that 
result from our simulations and those computed based on 
the EC model for the same wind magnetization. The value 
of Vts needed to compute the appropriate EC solution is 
derived again from Eq. E3 Actually, our simulations give 
a generally lower value of Vt s . Moreover this decreases as 



the ram pressure drops, eventually becoming negative (the 
shock starts collapsing back to the pulsar in the case with 
d = 0.003, n = 3 in Fig. [3). However, the comparisons in 
Figs. 01 and [5] arc not improved by using the exact values 
of Vt s given by the simulations. 

Again we find that the self-similar model gives a rea- 
sonably good approximation of the simulation results in 
the post shock region up to the radius at which the mag- 
netic pressure starts dominating over the thermal pres- 
sure, but fails, as could be expected, in the outer layers 
of the nebula. Moreover the singularity of the self-similar 
model, identified as the outer radius of the nebula to use 
for comparison with EC, is well inside the external bound- 
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Density 




R/R ts 

Thermal pressure 




R/Rts 

Magnetic pressure 




Fig. 6. Comparison of the results in the n = 2 case at 
t=1500 yr (dashed line) and in the case n = 3 case at 
t=760 yr (dotted line). The quantity p u that appears in 
the label of the y-axis is the density upstream of the ter- 
mination shock. Radial distances are normalized to the 
termination shock radius. In the different panels we show, 
from top to bottom, the density, thermal pressure and 
magnetic pressure profiles. In all cases a = 0.003. The 
solid line represents the EC model. 7 = 100. 



In Fig.|H|we compare the results obtained for a — 0.003 
and n — 2 and n = 3 at the time when the pulsar luminos- 
ity is the same. Despite the different conditions, the var- 
ious profiles coincide almost completely before the maxi- 
mum of the magnetic pressure, and deviations appear only 
in the outer layer of the nebula. This suggests that the in- 
ternal structure of a PWN is much more affected by the 
instantaneous properties of the wind, rather than by the 
overall pulsar history. Therefore, the appearance of the 
inner region of the PWN could in principle be used to es- 
timate the value of a from non-thermal emission, without 
sensitivity to the spin-down process. Once the magneti- 
zation of the wind is known, the ratio R c d/Rts could be 
used to infer the effect of spin-down (and eventually to 
determine n and r). 

4.3. Spin-down factorization 

In Fig. [7| the evolution of the ratio R c d/Rts is shown for 
the various values of a and n employed. Neglecting the 
spin-down effect may lead to considerable errors in the 
description of the PWN, even for ages comparable or less 
than the characteristic spin-down time. 

We want to stress that the ratio R c d/Rts plays an 
important role, being commonly used to infer the wind 
magnetization and from this the PWN properties. If the 
effects of the spin-down are not taken into account, the 
magnetization of the wind is overestimated, and the pro- 
files derived for the dynamically important quantities may 
be wrong. 

If the value of the magnetization parameter a is not 
known from independent constraints than the radius of 
the termination shock, to apply consistently the analytic 
model described in Section [21 the values of R c d and Rt s 
must be known at two different times, in order to fix both 
the normalization and a (see end of Section |2J). However 
in Subsection 14 . 21 we have verified that the EC model pro- 
vides a good description for R c d/R ts in the case n = 
or equivalently (<t. One can then use the result of the 
EC model to make up for the lack of observations over 
long time lapses. Thus the problem of evaluating a and 
the termination shock evolution could be reduced to an 
eigenvalue problem to be solved iteratively. 

From our simulations (see Fig.EJ) we find that the best 
fit for the ratio R c d/Rts is obtained multiplying the value 
given by EC (with Vt s given in Eq. 122(1 by 



ary as determined from the simulations. The outer part of 
the nebula shows also a positive velocity gradient. This 
is the effect of the extra energy in the outer layer: this 
layer was created when the pulsar was more energetic and 
carries more energy than it would if the pulsar luminos- 
ity had stayed constant, so that it tends to expand and 
causes the material more recently injected by the wind to 
be confined at smaller radii. The radius of the EC singu- 
larity corresponds, within a 10 % uncertainty, to the point 
where the speed of the flow is at its minimum. 



Co 



(l+f/r)"-(l + t/r) 
(n - l)i/r 



0.17+33 a 



(24) 



C a compensates for small differences that are present even 
in the case n = 0, and its values are: C Q = 1.04 in the HD 
case, and 1.06 and 1.07 for a = 0.0016 and a = 0.003 
respectively. Notice that the effects of the magnetization 
are summarized in the exponent. It is possible in principle 
to use this equation together with Eq.[53to estimate either 
a if n and r are known, or the spin-down parameters if 
the magnetization is independently known (for example 



10 



N. Bucciantini et al.: The effects of spin-down on the structure and evoiution of PWNe 



200 400 600 800 1000 1200 1400 



u = 0.0016 n = 




200 400 600 800 1000 1200 1400 



= 0.003 n = 2 




Fig. 7. Evolution of the PWN size, in units of the termi- 
nation shock radius. The upper three panels refer to the 
HD case, the middle three to a = 0.0016, and the lower 
three to a = 0.003. In each panel, the points refer to the 
result of our simulation; the dotted line is obtained using 
the value given by EC for the radius at which v = V c d\ 
the dashed line is the total size of the nebula in the EC 
model; the dot-dashed line is the fit done according to the 
correction in Eq. [31 the solid line, finally, represents the 
solution of Eq. ^3 normalized to the numerical results at 
time t — 50 yr. 

from the comparison of synchrotron and Inverse Compton 
emission) . 



Neglecting adiabatic losses, the formula obtained for 
a generic n is in shape similar to Eq. 1241 but with an 
exponent 0.5. Moreover Ea. 1241 will be valid only for small 
values of a given the fact that the exponent cannot exceed 
the value 0.5 (no adiabatic losses). More generally as a 
increases the value of the exponent will tend to 0.5. 

4.4. The Crab Nebula 

As a special case we shall consider the Crab Nebula. This 
is surely the best studied PWN: a wealth of information is 
available but still no answer has been given to a number 
of problems. 

We want to point out that ID models arc an over- 
simplification in general. More so for the Crab Nebula, 
where very clear axisymmetry is observed in the X-rays 
(the renowned jet-torus structure). This has suggested 
that important deviations from the spherically symmet- 
ric approximation might occur in the pulsar wind region, 
leading to the formation of a turbulent flow in the nebula 
(e.g. [Komissa rov fe Lyubarsky| 2003) . 

In modeling the Crab Nebula, one important pa- 
rameter is the actual value of the contact dis- 
continuity velocity; recent estimates IjFesen et al. 19971 
ISankrit fe Hester 19971 ISollerman et al. 2000jl give the 
value 1500 km/s with an upper limit of 2000 km/s. If 
we use the EC model without spin-down with the pre- 
scription in Ea. 1221 we find that, in order to end up at the 
present time with a nebula that has R c d/Rts = 20, a wind 
with a = 0.0009 or a — 0.0015 is required, depending 
on whether V c d — 1500km/s or V c d — 2000km/s is used. 
To include the effects of the spin-down, we use the fol- 
lowing values: r = 720 yr, n = 2.318 HCamilo et al. 20001 
|Lyne et al. 1993| |Lyne fc Manchester 1988|). U sing these 
values and the correction form given in Eq. 1241 we find for 
the magnetization of the Crab wind the value a — 0.0004, 
or, if the upper value of the expansion speed is used, 
a = 0.0009. 

We emphasize that these should be taken as estimates 
of the magnetization of the equatorial part of the outflow 
from the Crab pulsar, rather than being interpreted as 
reliable estimates of the overall magnetization of the Crab 
pulsar wind. 

5. Conclusions 

Our results show that, in the spherically symmetric ap- 
proximation, the EC model gives a good description of the 
internal profiles of PWNe when the effects of spin-down 
can be neglected (i.e. when the central object is young and 
not very powerful). Noticeable discrepancies arise only at 
the outer boundary of the nebula, where the initial con- 
ditions and the evolutionary history of the system both 
play a major role. However, the region just behind the 
termination shock and further out, up to the position of 
the equipartition point, is well reproduced. This is the re- 
gion of the nebula from which the high energy emission 
originates. 
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Even if the radial profiles of the various quantities are 
in good agreement, the ratio R c d/Rts inferred from the EC 
model by matching the boundary velocity for a given a is 
generally underestimated, due to the discrepancies at the 
boundary. Hence, if the ratio R c d/R ts estimated in this 
way is compared with the observed value to evaluate the 
wind magnetization, one generally overestimates the value 
of a. A better agreement is obtained if R c d is taken equal 
to the radius at which the singularity of the EC model is 
found. 

There is a good agreement in the HD case while a small 
difference, less than 10%, is found in the magnetized cases, 
probably because in these cases V c d ~ V asy . The EC model 
can be adapted to the evolution of V c d since, with respect 
to the sound speed crossing time in the PWN, variations at 
the border are very slow: changes of V c d are compensated 
by analogous variations of V ts . 

More generally, we expect that deviations from the 
spherically symmetric approximation might occur in the 
pulsar wind region, leading to the formation of a turbu- 
lent flow in the nebula that only multidimensional simula- 
tions can handle. In spite of being a simplification of the 
problem, ID models can help understanding the impor- 
tance of various processes like spin-down, mass loading or 
synchrotron losses, and to clarify how to keep them into 
account when comparing models with observations. 

Once the spin-down is included, the ratio R c d/Rts 
increases with respect to the constant luminosity case. 
Again, if this effect is not taken into account, the pulsar 
wind magnetization inferred from the PWN size might 
well be overestimated. With a simplified analytic model 
we are able to reasonably reproduce the ratio R c d/Rts 
given by the simulations. We also found a fit to the re- 
sults of our simulations in the form of a correction coeffi- 
cient that should multiply the standard EC expectation. 
However, the radial profiles of the various quantities given 
by the EC model fail to reproduce the structure inside the 
PWN beyond the point where magnetic pressure reaches 
its maximum. 

Given our fit, if the wind magnetization is known from 
the analysis of the non-thermal emission of the PWN, the 
size of the nebula can be used, in principle, to estimate 
n and r. Viceversa, if the pulsar spin-down properties are 
known, from the same expression it is possible to estimate 
a. 

As shown in the case of the Crab Nebula, inclusion 
of the spin-down effect reduces the estimated a by about 
one order of magnitude with respect to the value found 
by KC. We want to stress that this is not just a dynam- 
ical problem related to the size of the nebula alone. A 
less magnetized wind creates a less magnetized nebula. 
Synchrotron losses are less efficient and high energy par- 
ticles have a longer life and can move farther from the 
termination shock. This suggests that the interpretation 
of synchrotron maps should take into account the new re- 
sults we have presented. 
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